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Abstract 

We study analytically the dynamics of a network of sparsely connected in- 
hibitory integrate-and-fire neurons in a regime where individual neurons emit 
spikes irregularly and at a low rate. In the limit when the number of neu- 
rons — > oo, the network exhibits a sharp transition between a stationary 
and an oscillatory global activity regime where neurons are weakly synchro- 
nized. The activity becomes oscillatory when the inhibitory feedback is strong 
enough. The period of the global oscillation is found to be mainly controlled 
by synaptic times, but depends also on the characteristics of the external in- 
put. In large but finite networks, the analysis shows that global oscillations of 
finite coherence time generically exist both above and below the critical inhi- 
bition threshold. Their characteristics are determined as functions of systems 
parameters, in these two different regimes. The results are found to be in good 
agreement with numerical simulations. 
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1 Introduction 

Oscillations are ubiquitous in neural systems and have been the focus of several 
recent studies (for reviews see e.g. Gray 1994, Singer and Gray 1995, Buzsaki and 
Chrobak 1995, Ritz and Sejnowski 1997). In particular, fast global oscillations in the 
gamma frequency range (> 30 Hz) have been reported in the visual cortex (Gray et al 
1989, Eckhorn et al 1993, Kreiter and Singer 1996), in the olfactory cortex (Laurent 
and Davidowitz 1994) and in the hippocampus (Bragin et al 1995). Even faster 
oscillations (200Hz) occur in the hippocampus of the rat (Buzsaki et al 1992, Ylinen 
et al 1995). In some experimental data, (see e.g. Eckhorn et al 1993, Csicsvari et al 
1998, Fisahn et al 1998) individual neuron recordings show irregular spike emission, 
at a rate which is low compared to the global oscillation frequency^. This raises the 
question of whether a network composed of neurons firing irregularly at low rates can 
exhibit fast collective oscillations, which theoretical analyses and modelling studies 
may help to answer. 

^email: brunel@lps.ens.fr 
^email: hakim@lps.ens.fr 

■^Laboratory associated with CNRS, Paris 6 and Paris 7 Universities 
Fast oscillations may be due in some cases to a synchronized subset of cells with high firing 
rates. The observation of cells with the required property has been recently reported in (Gray and 
McCormick 1996). 
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Previous studies of networks of spiking neurons have mostly analyzed, or sim- 
ulated, synchronized oscillations in regimes in which neurons behave themselves as 
oscillators, with interspike intervals strongly peaked around their average value (see 
e.g. Mirollo and Strogatz 1990, Abbott and van Vreeswijk 1993, van Vreeswijk ct al 
1994, Gerstner 1995, Hansel et al 1995, Gerstner et al 1996, Wang and Buzsaki 1996, 
Traub et al 1996). Several oscillatory regimes have been found with either full or 
partial synchronization. A regime particular to globally coupled systems has been 
described where the network breaks into a few fully synchronized clusters (Golomb 
and Rinzel 1994, van Vreeswijk 1996). In some simulations of networks with de- 
tailed biophysical characteristics, cells fire sparsely and irregularly during a global 
oscillation (Traub et al 1989, Kopell and LeMasson 1994, Wang et al 1995), but the 
complexity of individual neurons in these models makes it difficult to clearly under- 
stand of the origin of the phenomenon. The possible appearance of fast oscillations 
in a network where all neurons fire irregularly with an average frequency which is 
much lower than the population frequency therefore remains an intriguing question. 
It is the focus of the present work. 

Recurrent inhibition plays an important role in the generation of synchronized 
oscillations as shown by in vivo (McLeod and Laurent 1996) and in vitro experiments 
(Whittington ct al 1995) in different systems. This has been confirmed by several 
modelling studies (van Vreeswijk et al 1994, Gerstner et al 1996, Wang and Buzsaki 
1996, Traub et al 1996). It has also been recently shown using simple models that 
networks in which inhibition balance excitation (Tsodyks and Sejnowski 1995, Amit 
and Brunei 1997a, van Vreeswijk and Sompohnsky 1996) are naturally composed of 
neurons with low and irregular firing. Simulations (Amit and Brunei 1997b) have 
shown that, in one such model composed of sparsely connected integrate- and- fire 
(IF) neurons, the highly irregular single neuron activity is accompanied by damped 
fast oscillations of the global activity. 

In order to study the coexistence of individual neurons with low firing rates and 
fast collective oscillations in its simplest setting, we analyze in the present paper 
a sparsely connected network entirely composed of identical inhibitory IF neurons. 
Our aim is to provide a clear understanding of this type of synchrony and to precisely 
determine : 

- i) under which conditions collective excitations of high frequencies arise in such 
networks 

- ii) what controls the different characteristics (amplitude, frequency, coherence 
time,...) of the global oscillation. 

Simulation results are presented first which shows that the essence of the phe- 
nomenon is present even in this simple system. Both the neurons firing rates and 
the auto-corrclation of the global activity are very similar to those reported in (Amit 
and Brunei 1997b). 

We begin by presenting simple arguments which give an estimation of the firing 
rate of individual neurons and the frequency of the global oscillation and which 
lead to think that the global oscillation only appears above a well-defined parameter 
threshold. 

In order to make the analysis more precise and complete, we then generalize the 
analytic approach of Amit and Brunei (1997a) which was restricted to the compu- 
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tation of firing rates in stationary states. The sparse random network connectivity 
leads the firing patterns of different neurons to be only weakly correlated. As a con- 
sequence, the network state can be described by the instantaneous distribution of 
membrane potentials of the neuronal population, together with the firing probability 
in this population. We obtain the coupled temporal evolution equations for these 
quantities, the time-independent solution of which coincides with the stationary so- 
lution of (Amit and Brunei 1997a). 

A linear stability analysis shows that this time-independent solution becomes 
unstable only when the strength of recurrent inhibition exceeds a critical level, in 
agreement with our simple arguments. When this critical level is reached, the sta- 
tionary solution becomes unstable and an oscillatory solution develops (via a Hopf 
bifurcation). The time scale of the period of the corresponding global oscillations is 
set by a synaptic time, independently of the firing rate of individual neurons, but 
the period precise value also depends on the characteristics of the external input. 

The analysis is then pushed to higher orders. We obtain a reduced evolution 
equation describing the network collective dynamics. The effects coming from the 
finite size of the network are also discussed. We show that having a large but finite 
number of neurons gives a small stochastic component to the collective evolution 
equation. As a result, it is shown that cross-correlations in a finite network present 
damped oscillations both above and below the critical inhibition level. Below the 
critical level, the noise controls the oscillation amplitude which decreases as the 
number of neurons is increased (at a fixed number of connections per neuron). Above 
the critical level, the main effect of the noise is to produce a phase diffusion of the 
global oscillation. An increase in the number of neurons results in an increase of 
the global oscillation coherence time and in a reduced damping in average cross- 
correlations. 

Finally, the effect of some of our simplifying assumptions is studied. We shortly 
discuss the effect of allowing variability in synaptic times and number of synaptic 
connections from neuron to neuron. We also consider the effect of introducing a more 
detailed description of postsynaptic currents into the model. The technical aspects 
of our computations are detailed in several appendices. 

2 Description of the network and simulations 

We analyse the dynamics of a network composed of N identical inhibitory single com- 
partment integrate-and-fire (IF) neurons. Each neuron receives C randomly chosen 
connections from other neurons in the network. It also receives Cext connections 
from excitatory neurons outside the network (see Fig. |lD. We consider a sparsely 
connected case with e = C/iV -C 1. 

Each neuron is simply described by its membrane potential. Let us suppose 
that neuron i receives an inhibitory (excitatory) connection from neuron j. When 
the presynaptic neuron j emits a spike at time t, the potential of the postsynaptic 
neuron i is decreased (increased) by J at time t + 6 and returns exponentially to 
the resting potential in a time r which represents the integration time constant of 
the membrane. In this simple model, the single time 6 is meant to represent the 
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Figure 1: Schematic diagram of the connect ionsSaJibfi-it^work of neurons; each 
neuron (indicated as an open disk) receives C inhibitory connections (indicated as 
black) from within the network and Cext excitatory connections (indicated as grey) 
from neurons outside the network. 
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Figure 2: Comparison of the synaptic response characteristics in our model and in a 
more realistic model. The top trace shows the presynaptic spike. The middle trace 
shows the corresponding postsynaptic current (PSC). The bottom trace shows the 
corresponding postsynaptic potential (PSP) for a neuron initially at resting poten- 
tial. Full lines: our model, in which the synaptic current is described by a delta 
function a time 5 after the presynaptic spike. Dashed lines: a more realistic synaptic 
response, in which the PSC is described by an a-function with latency (transmission 
delay) tl and synaptic time constant ts {t — TLj exp{—(t — tl) /ts) / ts ■ Our synaptic 
characteristic time 6 can roughly be identified with the sum of latency and synaptic 



decay time, tl + ts- See the discussion in Section O 
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transmission delays but also and most importantly, the longer time needed to ob- 
tain the full hyperpolarization of the post-synaptic neuron corresponding to a given 
presynaptic spike. Therefore, finding the correspondence between 6 and the different 
synaptic time scales of a more realistic description needs some care. As pictorially 
shown in Fig. ^ S should roughly be identified to the characteristic duration of the 
synaptic currents. In the following, we thus refer to S, which plays a crucial role in 
the generation of global oscillations, as the "synaptic time". The correspondence be- 
tween 6 and the different synaptic time scales of a more realistic description is further 
elaborated in Section [4.3| where synaptic current of finite duration are considered. 

Mathematically, the depolarization Vi{t) of neuron i {i = 1, . . . ,N) at its soma 
obeys the equation, 

TV,it) = -V,it) + Rl.it) (1) 

where Ii{t) are the synaptic currents arriving at the soma. These synaptic currents 
are the sum of the contributions of spikes arriving at different synapses (both local 
and external). These spike contributions are modelled as delta functions in our basic 
IF model: 

RU{t) = rJ2J^,J:^it-tJ-^) (2) 

j k 

where the first sum on the r.h.s is a sum on different synapses (j — 1, . . . , C -|- Cext): 
with postsynaptic potential (PSP) amplitude (or efficacy) Jij, while the second sum 
represents a sum on different spikes arriving at synapse j, at time t = tj+6, where 
is the emission time of k-th spike at neuron j. For simplicity, we take PSP amplitudes 
equal at each synapse, i.e. J^j = J^xt > for excitatory synapses and J^j = —J for 
inhibitory ones. External synapses are activated by independent Poisson processes 
with rate Uext- 

A firing threshold 6, completes the description of the IF neuron : when Vi{t) 
reaches 6, an action potential is emitted by neuron i, and the depolarization is reset 
to Vr < 6 after a refractory period r^p during which the potential is insensitive 
to stimulation. A typical value would be r^p ~ 2ms. We are interested here in 
network states in which the frequency is much lower than the corresponding maximal 
frequency l/r^p ~ 500Hz. In this regime, we have checked that the exact value of 
Trp does not play any role. Thus in the following we set r^p to zero, for the sake of 
simplicity. 

The outcome of a typical simulation is shown in Figs. ^ Neurons are driven 
by the random external excitatory input above threshold; however, since feedback 
interactions are inhibitory, the global activity stays at rather low levels (about 5Hz 
for the parameters indicated in Fig. §). For weak external noise levels [aext = 
ImV), the global activity (total number of firing neurons in 0.4ms bins) is strongly 
oscillatory with a period of about 7 ms, as testified by Fig. ^C. On the other hand, 
increasing the external noise level strongly damps and decreases the amplitude of 
the global oscillation. Note that the global activity should roughly correspond to 
the local field potential (LFP) often recorded in neurophysiological experiments. On 
the other hand, even when the global activity is strongly oscillatory, individual firing 
is extremely irregular as shown in the rasterfile of 50 neurons. Fig. (above the 
LFP), and in the inter-spike interval histogram (to the right of the spike rasters). In 
each oscillatory event only a small fraction of the neurons fire. 
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Figure 3: Left: Time evolution of the global activity (LFP) during a 100ms interval 
of the dynamics of a network of 5,000 neurons (total number of firing neurons in 
0.4ms bins), together with spike rasters of 50 neurons, for different values of the 
external noise: aext — 5mV (A), 2.5mV (B), and 1 mV (C). Right: autocorrelation 
of the global activity (AC) and inter-spike interval (ISI) histogram averaged over 
1000 neurons, corresponding to the left pictures. Note the different time scales of 
AC and ISI in abscissa. Parameters: 6 = 20mV, Vr = lOmV, r = 20ms, S = 2ms, 
C = 1000, J = O.lmV, i^ext = 25mV. 6 



This oscillatory collective behavior is also shown by fast oscillations in the tem- 
poral autocorrelation (AC) of the global activity which are damped on a longer time 
scale (Fig. to the right of the LFP). It is also reflected in the cross-correlations 
(CC) between the spike trains of a pair of neurons, which are typically equal to the 
AC of the global activity. 

These simulation results raise several questions on the origin and characteristics 
of the observed oscillations. What is the mechanism of the fast oscillation? In which 
parameter region is the network oscillating? What are the network parameters which 
control the amplitude and the different time scales (frequency, damping time con- 
stant) of the global oscillation? How do they scale with the network size? The model 
is simple enough and an analytical study gives precise answers to these questions as 
shown in the following sections. 

3 An analysis of the network dynamics 

Several features simplify the analysis as noted in a previous study (Amit and Brunei 
1997a) of the neuron mean firing rates. First, as a consequence of the network sparse 
random connectivity (C -C A^), two neurons share a small number of common inputs 
and pair correlations can be neglected in the limit C /N — > 0. Second, we consider a 
regime where individual neurons have a firing rate v low compared to their inverse 
integration time 1/r and receive a large number of inputs per integration time r, each 
input making a contribution small compared to the firing threshold (J ^ ^)0. In this 
situation, the synaptic current of a neuron can be approximated by an average part 
plus a fluctuating gaussian part, and the spike trains of all neurons in the network 
can be self consistently described by Poisson processes with a common instantaneous 
firing rate v{t) but otherwise uncorrelated from neuron to neuron (that is, between 
t and t + dt, a spike emission has a probability v{t)dt of occurring for each neuron 
but these events occur statistically independently in different neurons) 

The synaptic current at the soma of a neuron (neuron i) can thus be written as, 

RIi{t)= fi{t) + a^r],{t) (3) 

The average part ii{t) is related to the firing rate at time t — 5 and is a sum of local 
and external inputs 

fi = fil + fiext with fil = -CJu(t - 6)t, fiext = C extJ extl'extT (4) 

Similarly the fluctuating part, a^/Tr|i{t), is given by the fluctuation in the sum of 
internal and external poissonian inputs of rate Cu and Cext^^ext- Its magnitude is 
given by 

= af + aj^t with ai = J\JCv{t - 5)r, aext = Jext\/CextJ^extT (5) 

^Typical numbers in cortex are C = 5000, r = 20ms, ly — 5Hz, J — O.lmV, 9 — 20mV so that 
Ci^T is typically several hundreds while 9/ J is of order 100 (Abeles 1991, Braitenberg and Shutz 
1991). In the simulation shown in Fig. | Ciyr - 100, 0/J 200. 
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and r]i{t) is a gaussian white noise uncorrelated from neuron to neuron, {f]i{t)) = 



Before describing our precise results, it may be useful to give simple estimates 
which show how the neuron firing rates, the collective oscillation frequency and the 
oscillatory threshold can be obtained from Eqs.(|]-|). 

Let us first consider the stationary case. The case of interest corresponds to 
fi < 6. When expression @ is used for the synaptic current, the dynamics of the 
neuron depolarization (|I|) is a stochastic motion in the harmonic potential {V — /i)^ 
truncated at the firing threshold V = 6. The neuron firing rate z/q is the escape rate 
from this potential. For a weak noise, it is given by the inverse of the time scale of 
the motion 1/r diminished by an Arrhenius activation factor. So, one obtains the 
simple estimate (up to an algebraic prefactor). 



This becomes a self-consistent equation for uq once and a are expressed in terms 
of z/q using Eq. (^^. The simple estimate (|^) is made precise below by following 
Kramers's classic treatment of the thermal escape over a potential barrier (Chan- 
drasekhar 1943). 

The origin of the collective oscillation can also be simply understood. An increase 
of activity in the network due to a fluctuation provokes an increase in the average 
feedback inhibitory input. Thus after a period of about one synaptic time the activity 
should decrease due to the increase of the inhibitory input. This decrease will itself 
provoke a decrease in the inhibitory input, and a corresponding increase in the 
activity after a new period equal to the synaptic time. This simple argument predicts 
a global oscillation period of about a couple of times the synaptic time 6, not too 
far from the period observed in the simulations. However, it does not seem to have 
been noted previously that a global oscillation of period 6 can in fact occur only 
if it is not masked by the intrinsic noise in the system. The resulting oscillation 
threshold can be simply estimated in the limit where S is short compared to the 
time scale of the depolarization dynamics. During a short time interval 6, a neuron 
membrane potential receives from the local network an average input of magnitude 
Cvq5J. The fluctuation in its membrane potential in the same time interval (due 



to intrinsic fluctuations in the total incoming current) is a^Jd/r. The change in the 
average local input can be detected only if it is larger than the intrinsic potential 
fluctuations. A global oscillation can therefore occur only when 



These simple estimations are confirmed by the analysis presented below and 
replaced by precise formulas. 

3.1 Dynamics of the distribution of neuron potentials 

When pair correlations are neglected, the system can be described by the distribution 
of the neuron depolarization PiV, t), i.e. the probability of finding the depolarization 



and {Vi{t)vj{t')) = 5^,j5{t-t'). 




(6) 
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of a randomly chosen neuron at V at time t. This distribution is the (normahzed) 
histogram of the depolarization of all neurons at time t in the large N limit N ^ oo. 
The stochastic equation (|l],§) for the dynamics of a neuron depolarization can be 
transformed into a Fokker-Planck equation describing the evolution of their proba- 
bility distribution (Chandrasekhar 1943) 

The two terms in the r.h.s. of (|^) correspond respectively to a diffusion term coming 
from the current fluctuations and a drift term coming from the average part of the 
synaptic input. a(t) and are related to uit — 6), the probability per unit time 
of spike emission at time t — 6, hj Eq. Note that the Fokker-Planck equation 

has been used previously in studies of globally coupled oscillators (Sakaguchi et al 
1988, Strogatz and Mirollo 1991, Abbott and van Vreeswijk 1993, Treves 1993). 

The resetting of the potential at the firing threshold {V = 9) imposes the ab- 
sorbing boundary condition P(6', t) = 0. Moreover, the probability current through 
6 gives the probability of spike emission at t, 

At the reset potential V = Vr, P(y,t) is continuous but the entering probability 
current imposes the following derivative discontinuity, 

^(V- t) - ^(V- t) = (9) 

At = — oo, P should tend sufficiently quickly toward zero to be integrable, i.e. 

lim P{y,t) = Q lim VP{y,t) = Q. (10) 

Last, P{y,t) is a probability distribution and should satisfy the normalization 
condition 

/ p{y,t)dv = i (11) 

J —oo 

3.2 Stationary states 

We first consider stationary solutions P(y,t) = PoiV). Time independent solutions 
of Eq. (0) satisfying the boundary conditions (||,y,^) are given by 



,,2 



P„(V0.2^exp^-^^jy^ e^«-^Je"-*. (12) 
with 
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Figure 4: The neuron firing rate vs aext'- simulation (o); solution of Eq. (|l^)(full line); 
solution of the approximate asymptotic form ( |15|) (dashed line). Others parameters 
are fixed as in Fig. 2 : r = 20ms, J = O.lmV, C = 1000, N = 5000, 9 = 20mV, 
Vr = lOmV, fiext = 25mV, 6 = 2ms. 



(in (|T2|), 9(x) denotes the Heaviside function, 6(x) = 1 for x > and 6(x) = 
otherwise). The normalization condition (|ll|) provides the self-consistent condition 
which determines 



with 
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(14) 



,yr 



Vr-fJ-O 



. In the regime {6 — ^q) ^ aQ, Eq. ( P^ ) becomes 



(g - /io) 

O'OV^ 



exp 



(15) 



In Fig. (p, the firing rates obtained by solving Eq. ([T^) and (|T3|) are compared with 
those obtained from simulations of the network. It shows an almost linear increase 
in the rates as a function of aext in the range 3-6Hz and a good agreement between 
Eq. (n) and the results of simulations. The asymptotic expression (^) is also rather 
close to the simulation results in this range of a. 



3.3 Linear stability of the stationary states 

We can now investigate in which parameter regime the time independent solution 
{PqCV),^^) is stable. To simplify the study of the Fokker-Planck equation (|^, it is 
convenient to rescale P, V and u by 



P 



2ri^o 



-Q, y 



V 



u = Uo{l + n{t)) 



(16) 
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y is the difference between the membrane potential and the average input in the sta- 
tionary state, in units of the average fluctuation of the input in the stationary state. 
n{t) corresponds to the relative variation of the instantaneous frequency around the 
stationary frequency. After these rescalings, Eq.(|^) becomes 

dQ ^d'Q d (dQ Hd'Q\ 

where G is the ratio between the mean local inhibitory inputs and ctq, and H is 
the ratio between the variance of the local inputs and the total variance (local plus 
external) : 

^ 7_.. .. t2_,, ^2 

(To do o-Q erg 

These parameters are a measure of the relative strength of the recurrent inhibitory 
interactions. 

Eq. (p^ holds on the two intervals —oo<y<yr and yr < y < ye- The boundary 



^ = 1 = ^ = 12 = T2"' ^^^) 



conditions on Q are imposed at ye = and yr = . Those on the derivatives 
of Q read, 

dQ + dQ l + n{t) 

^{ye.t) = -{y-,t) - -iy^,t) = - , ^ ^^^^ _ (19) 

The linear stability of the stationary solution is studied in detail in Appendix 
1^. This can be done in a standard way (Hirsch and Smale, 1974) by expanding 
Q = Qo + Qi + ■ ■ ■ and n = rii + . . . around the steady state solution. The linear 
equation obtained at first order has solutions which are exponential in time, Qi = 
exp(wt/r)Qi, rii ~ exp(w/r)ni, where w is a solution of the eigenvalue equation 
( |66D of the Appendix. The stationary solution becomes unstable when the real part 
of w becomes positive. 

When the synaptic time 6 becomes much smaller than r, the roots w of this 
equation become large. We consider the regime 5/r ^ 1 but 6/r ^ 1/C, which is 
the relevant case in simulations and correspond to the realistic regime. 6/t ^ 1/C 
is needed because otherwise the equations giving G and H become inconsistent with 
the condition tuq -C 1. At the oscillatory instability onset, w is purely imaginary w = 
iuc, where uJc/t is the frequency of the oscillation which develops. The eigenvalue 
equation takes in the limit 6/t 0, u ^ oo the form 

-l) + H] exp{-iuJc5/T) = 1. (20) 



In this limit, the instability line in the parameter space (G, H) is obtained paramet- 
rically as 

G = 




n = sm + cos 
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Figure 5: Left: instability line in the plane (if, G\J 5/t). Full line: instability line for 
parameters of Fig.|^, and 5 = O.lr. Long-dashed line: 6 = 0.05r. Short-dashed line: 
asymptotic limit 6/t ^ 0. The stationary state (SS) is unstable to the right of the 
instability line, where an oscillatory instability develops (OS). Right: instability line 
in the plane {fJ.ext,<^ext)- Full line: parameters of Fig. 2, and 6 = O.lr. Short-dashed 
line is constructed taking the asymptotic instability line in the plane {H, G^JS/t), 
and calculating the corresponding instability line in {next,o'ext) with 6 = O.lr. The 
stationary state (SS) becomes unstable above the instability line. The long dashed 
line shows the average (fiext) and the fluctuations {(Text) of the external inputs when 
the frequency of a Poissonian external input through synapses of strength Jex* = 
O.lmV is varied. For low external frequencies the network is in its stationary state. 
When the external frequency increases the network goes to its oscillatory state (OS). 



H is by definition constrained to be between and 1 (it is the ratio between local and 
total variances): H = corresponds to the limit of very large external fluctuations, 
(Text ^ CT/, while H = 1 corresponds to cfext = 0. We find that the frequency of the 
oscillation varies from 

— = ^ when H = 0, to 
r 45 

'± = — when H = l. (21) 
r 2o 

This corresponds to an oscillation with a period between 85/3 and 46, not too far 
from the value 26 obtained by simple arguments. At the same time the critical value 
of G goes from 




when H = 0, to 
when H = 1. 



Again we find that it is proportional to yT/6 as anticipated. 

This instability line can be translated in terms of the parameters fiext, <^ext-, and 
calculated numerically using Eq. ( |5BD for any value of the network parameters. This 
line of instability in the plane {Hexti Cext) is shown in the right part of Fig. ^ The 
stationary solution is unstable above the full line. Thus, if the external input is 
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Poissonian, an increase in the frequency of external stimulation will typically bring 
the network from the stationary to the oscillatory regime, as indicated by the dashed 
line in Fig. ^ which represents the average (fiext) and the fluctuations (aext) of the 
external inputs when the frequency of a Poissonian external input through synapses 
of strength J^xt = O.lmV is varied. 



3.4 Weakly non-linear analysis 

The linear stability analysis of the previous section shows that a small oscillation 
grows when one crosses the instability line in the plane fiext, ^^ext- But it does not 
say much on the resulting characteristics of the resulting finite amplitude oscillation. 
In order to describe it and to be able to quantitatively compare analytic results to 
simulation data, one needs to compute the non linear terms which saturate the in- 
stability growth. This can be done in a standard manner (Bender and Orszag, 1987) 
by computing terms beyond the linear order in an expansion around the stationary 
state. The explicit computation is detailed in Appendix p. The collective oscillation 
is determined by the deviation of the neuron firing rate from its stationary value: 

= hxit) exp(zco'c^/'^) + exp(— zcjct/r) 

hx determines the amplitude of the collective oscillation as well as the nonlinear 
contribution to its frequency in the vicinity of the instability line. 

The analysis shows that the dynamics of the (small) deviation around the sta- 
tionary firing rate can be described by the reduced equation 

dn^ , ,o ^ s 

T—^ = Anx-B\hxYnx (22) 

in which A and B are complex numbers. The value of A comes from the linear 
stability analysis. If Re (A) < a small initial value of n\ decays and the stationary 
state is stable. On the contrary, if Re{A) > a global oscillation develops. When 
|ni| grows, the second nonlinear term on the r.h.s. of ( ^21) becomes important. It is 
found here that Re(i?) > (a "normal" or "supercritical" Hopf bifurcation) so that 
the nonlinear term saturates the linear growth. The characteristics of the oscillatory 
final state comes from the balance between the two terms. 

The explicit expression of A and B is given in Eqs. ([91|j9^) as a ratio of hyper- 
geometric functions of the network parameters. A depends linearly on the deviation 
of the parameters G and H from their critical values, i.e. G — Gc, H — He- In the 
limit 6/t ^ 0, the expressions of A and B simplify. For example, when H = (large 
external fluctuations), we find in the limit 6/t ^ 



A = l(i±^^.I,l.35 + 0.29P^-« 



c 



B 




4 + 97r 



2 



Gc S Gc 
13-5^2 9-5^2 /13-5V2 9-5^2' 

h I 1 

10 IStt V IStt 10 , 



(0.53 + 0.30i) (23) 
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Generally, the complex numbers A and B can be written in terms of their real 
and imaginary parts, A = Ar + iAi, B = Br + iBi. On the critical line, i.e. for 
G = Gc, H = He, Ar = Ai = 0; above the critical line an instability develops, 
Ar > proportionnally to G — Gc and H — Hf.. The amplitude of this instability is 
controlled by the cubic term. The stable limit cycle solution of Eq. (^), above the 
critical line, is 

hi{t) = Rexp (^iAuj-^ (24) 

where 

it='i/-— and Au = Ai — Bi — 

y JDr IJr 

The autocorrelation (AC) of the global activity, normalized by z/q, is, when Ar > 

0, 

C{s) = \im r \l + ni{t)){l + ni{t + s))dt (25) 

T^oo I — S JO 

= 1 + 2R'^ COS [{ujc + Auj)s/t] 

The AC is a cosine function of frequency {ujc + Auj)/t and amplitude R^. Compared 
with the AC function observed in the simulation. Fig. ^C, we see a qualitative dif- 
ference: there is no damping of the oscillation. The next Section shows that the 
damping is due to finite size effects. We analyze them before comparing quantita- 
tively the analytical results with simulations. 



3.5 Finite size effects and phase diffusion of the collective 
oscillation 

We discuss the effect of having a large but only finite number of neurons in the 
network. It is well-known that for stochastic dynamics, a sharp transition can only 
occur in the limit N ^ oo and that it will be smoothened by finite size effects. In 
the sparse connectivity limit, which allows to treat the quenched random geometry 
of the lattice in an annealed fashion^ the fluctuations in the input of a given neuron 
i can be seen as the result of the randomness of two different processes: the first is 
the spike emission process S{t) of the whole network; and the second, for each spike 
emitted by the network, is the presence or absence of a synapse between the neuron 
that emitted the spike and the considered neuron: if a spike is emitted at time t, 
Pi{t) = 1 with probability C/N, and otherwise. The input to the network is then 

Rliit) = -JTpi{t)S{t-5) 

Both processes can be decomposed between their mean and their fluctuation, 

p,(t) = ^ + 5p,(t), S{t) = Nu{t) + 6S{t) 

^Here we do not consider the correlations due to the quenched connectivity for finite e. These 
correlations would give small corrections to the parameters calculated in the limit e ^0, but do not 
give rise to qualitatively new effects for the global activity such as the phase diffusion phenomenon 
discussed in this section. 
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Thus the input becomes 

RIi{t) = fx{t) - jTNu{t)Spi{t) - JT^SS{t) 

in which is given by Eq. (^. The input is the sum of a constant part fi, 
and of two distinct random processes superimposed on /x: the first is uncorrelated 
from neuron to neuron, and we have aheady seen in Section ^ that it can be de- 
scribed by uncorrelated Gaussian white noises ay/rrjiit), i = 1,...,N where 
< rii(t)rij{t') >= 6ij6(t — t'). The second part is independent of i: it comes from the 
intrinsic fluctuations in the spike train of the whole network which are seen by all 
neurons. This part becomes negligible when e = C/N ^ 0, but can play a role as we 
will see when C/N is finite. The global activity in the network is essentially a Pois- 
son process with instantaneous frequency Nv{t). Such a Poisson process has mean 
Nv{t), which is taken into account in /i, and variance Nv{t)5{t — t'). The fluctuating 
part of this process is well approximated by a Gaussian white noise y/NuoC, (t) , where 
^{t) satisfies < ^{t) >= 0, < ^(t)^(t') >= 6{t - t'). Note that for simplicity we take 
the variance of this noise to be independent of time, which is the case for ni{t) -C 1. 
These fluctuations are global and perceived by all neurons in the network. Thus, the 
mean synaptic input received by the neurons becomes 



CJruit) + J^JeCu^T^iit) + 

l-^ext 

Inserting this mean synaptic input in the drift term of the Fokker-Planck equation, 
we can rewrite Eq. (0) as 

= l^f^ + - + ^v^e(t)]Q} + \^ (26) 

where r] denotes the intensity of the noise stemming from these global fluctuations. 
7] tends to zero as the network size increases 

V = Ve^ (27) 

Taking into account this global noise term in the derivation of the reduced equa- 
tion, we obtain, after some calculations described in Appendix |C|, 

= Ani- B\ni\'^ni + D^/^C{t) (28) 

in which A, B and D are given by Eqs. (|91| , |92|j9^) , and C is a complex white noise 
such that < C(t)C*{t') >= 5{t—t'). D is proportional to r], i.e. to both the square root 
of the connection probability and to the ratio between local and total fluctuations. 

Thus, the effect of the finite size of the network is to add a small stochastic 
component to the evolution equation of ni, Eq. (P^. Its main effect is to produce a 
phase diffusion of the collective oscillation []his global phase diffusion in a network of 
finite size is well-known (see e.g. (Rappel and Karma, 1996) for a simple example) 
which leads to the damping of the oscillation in the autocorrelation function. 
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Amplitude of the autocorrelation 



From the reduced Eq. (|2^), one can compute exactly the autocorrelation at zero time 
C(0) as shown in Appendix y. This gives : 

• In the stationary regime far from the critical line, Ar < 0, <^ 1: 

^(0) -^-^-^-O(^) (29) 



\Ar\ \N 

The amplitude of the fluctuations in the global activity are proportional to 
C/N and thus vanish when the connection probability goes to zero. 

• On the critical line, Aj. = 

The amplitude of the fluctuations are proportional to the square root of the 
connection probability. 

• In the oscillatory regime far from the critical line, Ar > 0,\D\/Ar <^ 1 : 

2 A 

C(0)-l~-^~O(l) (31) 

In this regime the amplitude of the oscillation is to leading order independent 
of the noise amplitude. 

Oscillations below the critical line 

In the stationary regime far from the critical line, the fluctuations of activity ni 
provoked by the noise term can be considered small and thus we can neglect the 
cubic term. It is then easy to calculate the autocorrelation (AC) of the activity, 

^^'^ = ^ + ^""""P (-^) ([^^ + ^^l^) (32) 

It is a damped cosine function. The damped oscillation has frequency {uJc+Ai) / t and 
damping time constant proportional to r/|Ar|. The amplitude of the autocorrelation 
function is proportional to C/N. 

Oscillations above the critical line 

In the oscillatory regime far from the critical line, we find in Appendix |C| an AC 
function of the form 

C{s) = 1 + 2^ cos ((c^e + Acu)s/r) exp ( -2!i£l ) (33) 
Br V 2 1 
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It is again a damped cosine function. The damping factor exp (— 7^(s)/2) is different 
from an exponential only at short times s ~ 5. At longer times, s ^ S, we obtain 
again an exponential 

The damping time constant is proportional to leading order in \D\ to 1/|-Dp ~ N/C, 
i.e. to the inverse of the connection probability. When goes to infinity at C fixed 
the 'coherence time' of the oscillation increases linearly with A^. 

This 'phase diffusion' effect is the main finite size effect above the critical line. 
Both the amplitude and frequency of the oscillation are essentially unaffected by 
these finite size effects. 



3.6 Comparison between simulations and theory 

The autocorrelation (AC) of the global activity was computed for each set of pa- 
rameters from a simulation of 20 seconds. Few longer simulations were performed 
as a check. The autocorrelation obtained in the longer simulations are essentially 
identical to the one obtained in the 20s simulation. 

Since the analysis predicts AC functions described by damped cosine functions, 
a least square fit of all AC functions was performed with such functions. Thus the 
full AC is reduced to three parameters, its amplitude at zero lag Co, its frequency 
a;, and its damping time constant (or coherence time) Tc 

C(s) = 1 + Co exp cos(c<Js) 

We then compared the result of the fitting procedure with the analytical expressions. 

We have varied the magnitude of the external noise aext from to 5 mV. This 
brings the network from the 'oscillatory' to the 'stationary' state. 

In Fig. ^ we plot together the results of simulations and theory. In these figures 
the diamonds are the simulation results; the dashed lines, the analytical results. In 
A, the short-dashed line indicates the amplitude in the limit N ^ oo, while the 
long-dashed line indicates the amplitude calculated analytically taking into account 
finite size effects. Last, the crosses are obtained simulating numerically the reduced 
equation, Eq. |2^. We find that, in the 'stationary' regime as well as in the oscillatory 
regime close to the bifurcation point, the amplitude of the oscillation obtained in 
the simulation is in very good agreement with the calculation (Fig. p.A). On the 
other hand, as the amplitude of the oscillation becomes of the same order as the 
average frequency, Co ~ 1, higher order effects become important and the calculation 
overestimates the amplitude of the AC. For the frequency of the oscillation (Fig. §.B), 
the calculation reproduces quite well the results of the simulations, except for very 
low noise levels, for which we are rather far from the bifurcation point. Note that 
the frequency ranges for this set of parameters from 70 to 180Hz, depending on the 
level of external noise. Thus, without varying the time constants r and 6, we find 
that the same network is able to sustain a collective oscillation at quite different 
frequencies. 



17 



A 




B 




c 



70 I r- 

60 h o 

--^ 
<> ^ 
50 - \ 



'rc(ms) 4Q ^ 
30 

20 h 
10 




1 1.5 2 2.5 3 3.5 4 4.5 5 



Figure 6: Parameters of the AC function vs (Text- A. Amplitude of the AC at zero 
lag. B. Frequency. C. Damping time constant. Diamonds: simulation of the full 
network. Crosses : simulation of the reduced equation. Dashed hues: theory. In 
A, the short-dashed line represents the amplitude in the limit N ^ oo Parameters: 
T = 20ms, J = -O.lmV, C = 1000, N = 5000, 9 = 20mV, K = lOmV, fx^xt = 25mV, 
5 — 2ms. 
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Last, the approximate analytical expressions for the damping time constant agree 
well with the simulation away from the bifurcation point, as expected (Fig. |^.C). On 
the other hand, the simulation of the reduced equation is in good agreement with 
the network simulations in the whole range of (Text- 

In Fig. 1^ we compare the full AC functions from theory (simulation of the reduced 
equation) and network simulations in three regimes, to show the good agreement 
between both. 

4 Extensions 

In the previous sections a very simple network has been analyzed and the question of 
the effect of some of our simplifying assumptions legitimately arises. In particular, 
we have chosen exactly identical neurons. It can be wondered how the results are 
modified when some variations in neuron properties are taken into account. In order 
to address this question, we show how the previous analysis can be generalized in two 
cases. Since we have seen that the oscillation frequency is tightly linked to synaptic 
times, the effect of a fluctuation in synaptic times is investigated first. We then 
consider the effect of a fluctuation in the number of connections per neuron which 
has been found to result in a wide spectrum of neuron steady discharge rates (Amit 
and Brunei, 1997b). In both cases, it is reassuring to find that the picture obtained 
from the simple model analysis remains accurate. We finally consider a model with 
synaptic currents of finite duration to analyse more precisely which time scale plays 
the role of our "synaptic time" in this more realistic case. 

4.1 Effect of inhomogeneous synaptic times 

The analysis can easily be extended to the case in which time constants at each 
synaptic site are drawn randomly and independently from an arbitrary probability 
density function (pdf) Pr(5) (see Appendix 0). In the following we consider the case 
of a uniform pdf between and 25. 

Fig. ^ shows how the instability line is modified by random synaptic times. The 
region where the oscillatory instability appears shrinks to the area above the dashed 
line. As the distribution of synaptic times widens, the stationary state becomes 
more stable. The introduction of random synaptic times also slightly reduces the 
frequency of the oscillation. 

The critical line is thus quite sensitive to the distribution of synaptic times. In 
fact, distributions of synaptic times can be found such that the stationary state is 
always stable (e.g. for an exponential distribution Pr((5) = exp(— 5/5o)/5)- 

4.2 Effect of inhomogeneous connectivity 

The analysis can also be extended to the case when the number of connections 
impinging on a neuron is no longer fixed at C, but rather connections are drawn at 
random independently at each site. In that case the number of connections received 
by a neuron is a random variable with mean C and standart deviation ~ \fC. 
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Figure 7: AC for: A. a^xt — 2mV. B. a^xt — 3mV. C. o'ext = 4mV. Parameters as 
in Fig. 1^. Full lines: network simulation. Dashed lines: theory (simulation of the 
reduced equation). 
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Figure 8: Instability line in the plane {/J^ext, <^ext) for r = 20ms, J = O.lmV, C — 1000, 
9 — 20mV, Vr — lOmV, 5 — 2ms. Full line: all synaptic times equal to 5. Dashed 
line: synaptic times drawn from a uniform distribution from to 25. 
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Figure 9: Effect of inhomogcncity in the connections on the instability line in the 
plane {i^exu (^ext) for r = 20ms, J = -O.lmV, C = 1000, 9 = 20mV, Vr = lOmV, 
S — 2ms. Full line: all neurons receive C connections. Dashed hne: connections are 
drawn randomly and independently at each synaptic site with probability C/N. 
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Figure 10: Left: Distribution of spike rates (Histogram: simulation. Dashed line: 
theory). The distribution is similar to a Gaussian, unlike the distributions observed 
in (Amit and Brunei 1997b), which are much wider, due to the balance between 
excitation and inhibition. Right: Relative amplitude of CC between individual neu- 
rons and the global activity vs neuronal firing rate (Diamonds: simulation. Full line: 
theory), r = 20ms, J = -O.lmV, C = 1000, 6 = 20mV, K = lOmV, 5 = 2ms, 
fiext = 25mV, aext = 2.58mV. 



This inhomogeneity in the connectivity provokes a significant inhomogeneity in the 
individual spike rates even for C large, because differences between the average 
input received by two neurons are of the same order as the SD of the synaptic input. 
The distribution of frequencies for an arbitrary network of excitatory and inhibitory 
neurons has been obtained in (Amit and Brunei 1997b). The main steps leading to 
this distribution are described in appendix 0. Next we study how inhomogeneity 
affects the dynamical properties of the network. Fig. ^ shows that the instability line 
is almost unaffected by the inhomogeneity. The frequency of the global oscillation 
is also very close to the one of the homogeneous case. 

Amit and Brunei (1997b) had shown by simulations that the degree of synchro- 
nization of a neuron with the global activity is strongly affected by its spike rate: 
neurons with low firing frequencies tend to be more synchronized with the global 
activity than neurons with high frequencies. In appendix |^ we calculate analytically 
the degree of synchronization of individual neurons as a function of their frequency. 
The result is shown in Fig. |10| in which the relative amplitude C(z/) of the cross- 
correlation between neurons firing at frequency u and the global activity obtained 
analytically is compared with the result of simulations. It shows indeed that low-rate 
neurons are more synchronized with the global activity than high-rate neurons. The 
relative amplitude of the cross-correlation between two neurons of frequency i>i and 
z/2 is given by the product of the two amplitudes, C {1^1)0 {1^2). Note that the hetero- 
geneity in rates and cross-correlations is not very pronounced here, because near the 
critical line the fluctuations in the external input dominate the local fluctuations, 
which tends to suppress this heterogeneity. In a network with both excitatory and 
inhibitory neurons with an external excitatory input of the same order than the in- 
ternal excitatory contribution, this heterogeneity is much more pronounced (Amit 
and Brunei 1997b). 
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4.3 Effect of more realistic synaptic responses 

Our analysis has been carried out for synaptic currents which are described by a 
delta pulse. One may wonder how the analysis generalizes for more realistic postsy- 
naptic currents. We consider a function f{t) describing the shape of the postsynaptic 
current when a spike is emitted at time t = (see e.g. Gerstner 1995 for a review of 
different types of synaptic responses). f(t) is chosen such as 

dtfit) = 1 

An example often used in modelling studies and shown in Fig. ^ is the a-function 
with a latency tl and a characteristic synaptic time ts'- 



fit) 



^^exp(- 



■TS 







for t > Tl 
otherwise. 



(34) 



The total synaptic current arriving at neuron i is now 



In the diffusion approximation the synaptic current becomes 

RIi{t)= ^^{t) + E,{t) 

in which the average part is given as a function of the frequency v and the synaptic 
response function / by 

liit) = fi,^t -CJ j dt'u{t')f{t - t')T. 

On the other hand, the fluctuating part Hj(t) can no longer be approximated by 
a pure white noise and exhibits temporal correlations at the scale of the width 
of the PSC function f{t). These temporal correlations in the currents complicate 
significantly the analysis, since the evolution of the distribution of the membrane 
potentials is no longer given by a simple one- dimensional Fokker-Planck equation. 
For the case of the a-function, we would need to solve the problem described by 
a three dimensional Fokker-Planck equation. Such an analysis is beyond the scope 
of the present paper. Here, we choose to ignore, as a first approximation, these 
temporal correlations. Thus we consider only the effect of the PSC function on the 
average synaptic currents. In this approximation, the effect of the PSC function 
becomes equivalent to that of a distribution of synaptic times in the delta pulse PSC 



case considered in section [4.1| . For example, in the limit in which ts and tl are small 
compared to the integration time constant, the equations for the bifurcation point 
are 
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Figure 11: Dependence of the frequency of the oscillation near the bifurcation thresh- 
old on the synaptic decay time constant ts, for = 2ms. Network parameters as 
in Fig. 1^. External inputs have iie^t = 25mV, (Text = 2mV. This point is near the 
bifurcation line in the whole range of r^. o: simulations. Full lines: frequency given 
by the approximate analysis, Eq. ^ for if = 1 (lower curve), and H = (upper 
curve). 



In the case tl = (zero latency) the equations simplify to 

G = 2v4J— u; (36) 

r 

H = l-^uj^ + 2^LU (37) 

In the case H = 1, the frequency of the oscillation near the bifurcation point is equal 
to l/(7rrs). Note that the dependence of the frequency on ts in the a function PSC 
case is similar to the dependence on 6 in the delta pulse PSC case, Eq. (^Tf). 

To check the validity of this approximation, we have performed numerical simula- 
tions with fixed latency tl = 2ms, varying the decay time constant of the inhibitory 
post synaptic currents (IPSC) rs- The results are shown in Fig. |Tl]. The approximate 



analysis predicts the frequency is in the region between the two full lines (correspond- 
ing to if = and H = 1). Simulation results deviate from the approximate analysis 
already at rather small values of r^, because of the effect of temporal correlations 
in the synaptic currents, which have the same scale as the period of the oscillation. 
Nonetheless the approximation gives a good qualitative picture of the dependence of 
the frequency on ts- 

Note that the frequencies obtained in this way can be directly compared to the 
data of (Whittington et al 1995, Traub et al 1996) since the decay time constant of 
the PSCs can be identified with their parameter tgaba- The frequencies obtained 
in the simulations are very close to the ones obtained in that study. For example, 
we obtain a frequency of about 40Hz when ts = 10ms, in agreement with the in 
vitro recordings and the simulations of the more complex model of (Whittington et 
al 1995, Traub et al 1996). However, one has to be careful with such a comparison, 
since in that in vitro study, interneurons seem to fire at population frequency. 
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5 Conclusion 



We have studied the existence of fast global oscillation in networks where individ- 
ual neurons show irregular spiking at a low rate. We have first shown that the 
phenomenon can be observed in a sparsely connected network composed of basic 
integrate and fire neurons. In this very simplified setting, the phenomenon has been 
precisely analyzed. At the simplest level, it differs from other modes of synchronisa- 
tion which lead to global oscillation in that recording at the individual neuron level 
shows a stochastic spike emission with nearly Poissonian interspike intervals and 
little indication of the collective behavior (see the ISI histograms in Fig. This 
oscillation regime has some similarity with that obtained in Wang et al (1995) where 
a hyperpolarization-activated cation current seems to play the role of our random 
external inputs in generating intermittent activity in the network. This type of weak 
synchronization has sometimes been rationalized as coming from filtering of external 
noise by recurrent inhibition (Traub et al 1989 and refs. therein). Our analysis leads 
to a somewhat different picture. 

We have found that, in the limit of an infinite network, the global oscillation is 
due to an oscillatory instabihty (a supercritical Hopf bifurcation) of the steady state. 
This instability occurs at a well defined threshold and arises from the competition 
between the recurrent inhibition which favors oscillations and the intrinsic noise in 
the system which tends to suppress it. 

We have found that the global oscillation period is controlled by the synaptic 
time. This appears to agree with previous experimental findings on slices of the 
rat hippocampus and with simulations results (Whittington et al 1995, Traub et al 
1996) where it is however assumed that neurons fire at population frequency, unlike 
those of our model. A similar decrease in population frequency when the GABA 
characteristic time is varied is also observed in a recent in vitro experiment in which 
neurons fire sparsely (Fisahn et al 1998). More work is necessary to clarify the 
relative roles of the different time constants (latency, IPSO rise time, IPSO decay 
time) that are commonly used to describe the synaptic response. 

The oscillation period also depends on the characteristics of the external input, 
and particularly on the magnitude of the external noise, as shown by Fig. ^. The 
initial rise in the frequency when one increases aext followed by a saturation at 
sufficiently large aext looks in fact similar to the dependence of the frequency on 
the amount of glutamate applied to hippocampal CAl region in vitro (Traub et al 
1996). Our network is in a stationary state when external inputs are low and switches 
to an oscillatory regime when the magnitude of the external inputs is increased. 
This phenomenon resembles the induction of a gamma rhythm in the hippocampal 
slice mediated by carbachol (Fisahn et al 1998), and the induction of faster 200Hz 
rhythms, believed to be provoked by a massive excitation of CAl cells through 
Schaeffer collaterals (see e.g. Buzsaki et al 1992). It is also interesting to note that 
a single network, with its internal parameters fixed, is able to sustain collective 
oscillations in different frequency ranges, when the characteristics of the external 
input are varied. 

In a finite network, the sharp transition is smoothened but the global oscillation 
has different characteristics above and below the critical threshold. Below thresh- 
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old, its amplitude decreases as the network size is increased. Above threshold, an 
increase in the neuron number does not greatly modify the oscillation amplitude but 
increases its coherence time. It has been shown that the whole picture of a Hopf 
bifurcation with a well-defined threshold remains accurate when some of our simpli- 
fying assumptions are relaxed. It would be interesting to extend this finding to more 
realistic descriptions. 

Our analysis also raises the important question of the synchronisation mode used 
in real neural systems. Do neocortical or hippocampal neurons behave as oscillators 
with a frequency equal to the population frequency, or irregularly with firing rates 
lower than the population frequency? In hippocampus, pyramidal cells seems clearly 
to be in a irregular, low rate regime, during in vivo gamma (Bragin et al 1995), 
in vivo 200Hz (Buzsaki et al 1992) and in vitro gamma oscillations (Fisahn et al 
1998). More recent experimental data indicates that interneurons also typically fire 
at a lower frequency than the population frequency during 200IIz oscillations in 
CAl (Csicsvari 1998). Further experimental work is needed in order to clarify this 
important issue. 

We have obtained a reduced description of the collective dynamics. The analysis 
can certainly be extended to more complicated networks, composed of neurons of 
different types or that are spatially extended. This reduced description will hopefully 
prove useful in clarifying the mechanisms of long range synchrony and in studying 
propagation phenomena (Delaney et al 1996, Prechtl et al 1997). 

Finally, and most importantly, the exact roles of fast oscillations remain, at 
present, unclear. Are they useful for putting in resonance different neuronal popu- 
lations as it has been suggested? Can they serve to build a fast detector with slowly 
firing neurons? Are they used as a clock mechanism? Or do they reflect the use- 
fulness of having a network where different neuronal populations flre in succession 
on a short time scale, to code spatial information in the temporal domain? Recent 
experiments (MacLeod and Laurent 1996, Stopfer et al 1997) make us hope that 
elucidating the real meaning of these collective oscillations, at least in some neural 
systems, is now an attainable goal. This is a question to which we hope to return in 
the future. 
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Appendix 

The details of our computations are given in the following. We have found it conve- 
nient to use the rescaled variables. 



P=^g, G = ^ = if=^^ = 4^, (38) 



(7o (7o (7o CTg (Jq 
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V - fJ^O 6 - llQ Vr- llQ , , . . 

y = , ye = , yr = , i' = + n-it)) (39) 

(To (To (To 

J and G are positive. 

Using Eqs. (|38|j39D the Fokker-Planck equation, Eq. (|^) becomes 

r- = /:[Q] + .(t-5)(^G— + (40) 

where the hnear operator C is defined as 

1 d^O d 

The equation is vahd on the two intervals —oo<y<yr and i/r < y < ye- 
The boundary conditions at i/r and ye become: at ye 

at yr 

(the square bracket denotes the discontinuity of the function at y namely, = 
lime^o{/(l/ + e) ~ fhj"^)})- Note the term in the r.h.s. of Eqs. (p|,^) are identical. 
Thus, when we study the Fokker-Planck equation at different orders, we will mention 
only the condition at ye- The condition at y,. can be obtained by replacing the value of 
the corresponding function at ye by the discontinuity of the function at y^- Moreover 
Q{y,t) should vanish sufficiently fast ai y = — oo to be integrable. 
The steady state solution obeys 

C[Qo] = (43) 

and 

It is given by 

Qo(y) = |'"PK|£f '"^Kl ^^'^ (45) 
^"^^^ \ exp{-y^)jy^duexp{u^) y < yr ^ ^ 



From (|43| , |4^) , one easily obtains the values of higher derivatives of Qo at y = ye 
and their discontinuities at y = yr, which will be used in the following, using the 
recurrence relation 



27 



A Linear stability 

The function Q can be expanded around the steady state solution Qoiy) as 

Q{y) = Qo{y) + Qi{y, t) + Q2{y,t) + --- 

n{t) = ni{t) + n2{t) + ■ ■ ■ (47) 
At first order, one obtains the hnear equation 

= cm +n,{t- 6) [G^ + j (48) 

together with the boundary conditions 

Q^iye, t) = 0, ^(ye) = -^i(t) + Hn,{t - 6) (49) 

and 

[^C-=0' = -n.it) + Hn,{t - 6) (50) 



Eigenmodes of (]48|) have a simple exponential behaviour in time 

Qi{y,t) = exp(At/r) ni{X)Qi{y, A), ni{t) = exp(At/r) hi{X) 
and obey an ordinary differential equation in y 

XQ.iy, A) = C[Q,U A) + e'^'^^ (g^ + (51) 
together with the boundary conditions 

Q^iye^t) = 0,^(2/e) = -I + H exp{-\5 /r), 

and similar conditions at yr- 

The general solution of Eq.(pID can be written as a linear superposition of two 
independent solutions 0i_2 of the homogeneous equation 1/20" + ycj)' + (1 — A)0 = 
plus a particular solution which can be obtained by differentiating Eq. (^) with 
respect to y, 

n (v \] = f ^) + l^tWMy, A) + Qliv, A) y>yr . . 

^ 1 ar(A)0i(2/,A)+/?r(A)02(2/,A) + g?(2/,A) y < y^ ^ ^ 

with 

1^1 + A dy +2(2 + A) dy^ ) ^^^^ 

Solutions of the homogeneous equation l/20" + ?/0' + (1 — A)0 = can be obtained 
by their series expansion around y = 0. They are found to be a linear combination 
of two functions. The first one can be chosen as 

My, A) = 1 + E(-i)"^ n + 

n=l fe=0 ^ 
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It coincides with the confluent hypergeometric function M[(l — A)/2, 1/2, — y^] (see 
e.g. Abramovicz and Stegun 1970) . A second independent solution can also be 
expressed in terms of the hypergeometric function M as 

2.M (l - ^. |, -y^) . 2y . ^ " ^) (55) 

The asymptotic behaviour of both functions can conveniently be obtained from the 
following integral representations valid for Re(A) < 1/2 

My,^) = r^n-xs / dte~' cos{2yVi)t- — 

r(— ) Jo 

l--,-,Vj = f(7TTyX dte-' sm{2yVi)t-— (56) 



(after replacing the cosine and sine in (|56|) by their series expansions it is easily 
checked that the obtained series in powers of coincide with (|5^) and ([55|)). The 
following asymptotic behaviours are found for y ^ —oo 



We find it convenient to choose 02(z/,i^) as the particular combination of these two 
functions which decays exponentially (i.e. like \y\~^ exp{—y'^)) a.t y = — oo, 



My, ^) = -#T M^-^X -A + 2!,M ( 1 - ^, 5 (59) 



r(i±i) \ 2 -i- " )• r(A 



2 2' 



Thus for Qi{y,t) to be integrable on [—oo,yg] we need to require ai = in (|52|). 
For further reference, we give the asymptotic behaviour for A2 = Im(A) +00, 



MyAi + i^2) ~ cosh[?/JA2-iAi(l + i)]exp(-yV2) (60) 



02(2/,Ai + zA2) ~ -^g^exp[yyA^^^(l + z)-yV2] (61) 

I 2 j 

where the determination of the square root is fixed by requiring it to be positive for 
Ai = 0. 

Finally, we note that the Wronskian Wr of 0i and 02 obeys the first order equation 
Wr' = —2yWi and has therefore the simple expression 

Wr(0i, 02) = 0102 - 0'i02 = exp(-y2) (62) 

(the prefactor being fixed by (^,|6T|)). 
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The four boundary conditions (HUlBUp give a linear system of four equations for 
the four remaining unknowns ai,ai,Pi and The condition a{ = needed 
to obtain an integrable Qi{y,t) gives the eigenfrequencies of the hnear equation 
(^81). To obtain the required solvabihty condition and the alhed solutions, we find it 
convenient to use first the two boundary conditions ( ^9]) to obtain ai and f3i. This 
gives 

1 



{Mye){l-He-''/^)-W2[Ql] (ye)} 



fit 



,{y,){l-He-^'l^)-WAQf\ 



(63) 
(64) 



Wiiye) 
1 

~WrM 

where Wr denotes the Wronskian of 0i and 02, Eq- (0), and Wj{j = 1,2) the 
Wronskian of the function in its argument and 0i,2 

[q] = Qcj)'. - Q'ct)j for 3 = 1, 2. 
For matters of convenience we define 0i^2 and Wi^2 by 



yi,2 
Wr ^ 



1^1,2 Q? 



Wr 



The two boundary conditions at y = y,. 
Pi — Pi with ye replaced by yr 



give similar equations for ai — and 



Pi 



at-Hyr){l-He-^"^) 



W2 


q{ 


(y) 


Wi 




iy)_ 



Vr 
Vr 
Vr 



(65) 



The two expressions (|63|j65|) together with = give the solvability condition and 
the equation for the eigenfrequencies of (HSf) 



Wo 



Ql {y 



Wo 



Ql {y) 



Vr 



(66) 



When the synaptic time 5 becomes much smaller than r, the roots A of this 
equation become large. Considering for definiteness roots A = Ai + 1X2 with A2 > 0, 
in the limit |A| — > +00, A2 +00, one obtains from (|61|) that dy(t>2{ye) ^ dy(j)2{yr) 
and dy(j)2{ye) ~ \/\2 — iXiil + i)02- We then note that for Eq. (|HB| ) to have such 
a root, we need G ~ v |A|. Since if < 1 by definition, we can neglect the terms 



(67) 



proportional to H in Q\ and finally obtain 



G- 



~XS/t 



A 



We focus on the root with the largest real part (together with its complex conjugate). 
Its real part becomes positive, A = zA2 = ioJc when 

;i - i)Ge~^^^^l^ 



I.e. 



1 - He-''^^^/^ + 



H 







G = ^JZ^c sin (tUc^ / r) 

sin {ujc5 /t) + cos {ujc5 /t) . 
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B Weakly non- linear analysis 



Our aim is to determine the lowest non-linear terms which saturate the instability 
which appears when one crosses the critical line in the plane ^ext, o'ext- This deter- 
mines the amplitude of the collective oscillation as well as the nonlinear contribution 
to its frequency in the vicinity of {Gc, He) . We follow the usual strategy of pushing 
the development ( P7| ) to higher order. One finds that the nth-order term obey in- 
homogeneous linear equations with forcing terms formed by quadratic combinations 
of lowest-order terms. We first determine the second-order terms which are forced 
by quadratic combination of first-order terms and therefore oscillate at and 2uJc- 
At third order, the coupling between first and second order term generate forcing 
terms at ujc and StUc- While there is no problem to determine the 3ujc contribution, 
the Uc forcing is resonant and generates secular terms. The dynamics of the first- 
order terms amplitude is determined by the requirement that it cancels the unwanted 
secular contribution. The computation is not specially difficult but rather long. 

We substitute the developments (^) of Q{y,t) and n(t) in Eq. (|T^ anticipating 
that the development parameter is of order of the square root of the differences 
G — Gc, H — He- Departure of G from Gc and of H from will therefore only affect 
the third-order terms. 

The first-order terms have already been obtained, 

Qi{y,t) = e'^^'/^niQi{y, iujc) +C.C. 

n^(t) = e*"=*/"r2i(icjc) +C.C. (68) 



where Qi is given by Eqs. (^,^,^, In Eq. (^), we recall that c.c. means 
that complex conjugate terms to those explicitly written have to be added. In the 
following, we omit the explicit mention of the variable A to lighten the notation since 
functions of A will all be evaluated at iuc (except when explicitly specified otherwise). 

By differentiation of Eq. (^8]), one can easily obtain recursively the values of 
higher derivatives of Qi at y = yg and their discontinuities at ?/ = yr-, which will be 
used in the following. 

B.l Second order 

We first determine the second-order terms. They obey the equation 

= + - S) [G.^ + j + n,(t - 5) [G.^ + 

together with the boundary conditions 

Q2{ye, t) = 0, ^{ye) = -n^it) + Hn^it -6)- H^n\{t -5) + Hn,{t)n,{t - 5) (70) 
dy 

and a similar condition in yr. 

From (|68D, the forcing term on the r.h.s of Eq. (^) contains terms at frequencies 
2uj(. and 0. Therefore, we search Q2{y,t) and n2{t) under the form 

Q2{y,t) = e''-^'/^nlQ2,2{y)+e-''-^'/^{nl'f^^^^^ (71) 

n2{t) = e^'^^"^hlp2,2 + e-^'^^"^{n\fp\2 + \ni\^P2fl (72) 
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Substitution of ( [72| ) into ( |5UD shows that (^2,2 obeys the ordinary differential 
equation 



dy 2 y 
together with the boundary conditions 



(73) 



Q2,2{ye,t) = 0, ^^^{ye 



dy 

and a similar condition in 

As above, the general solution of ( |75D is written as a superposition of solution of 
the homogeneous equation and a particular solution 



Q2Ay) 



atMy^ 2iitjj + My, 2ituc) + P2,292?2 + 92^2 y>yr 
(^2 My, 2it^c) + 1^2 My, 2«c^c) + p2,2Qf2 + Qt^ y <yr 



(74) 



where 



/\so ~2iuJcS/r 

V2,2 ~ 



l + 2iujc dy A{1 + iuc) dy"^ 



Q2 2 can be obtained by differentiation of Qo and Qi using (^) and (|51| ) and involves 
only terms of lower order which have already been determined, 



QUy) 



-2iLUci/T 



\l + iuJc dy 2(2 + iujc) dy'^ J 

d'^Qo HcGc d^Qi] 



Hi 



d'Qo' 



^2(l + zcJc)2 t/^2 2{1 + iuJc) {2 + iuJc) dy^ 8(2 + zcjj^ dy 



The four boundary conditions for Q2 determine the four unknowns at , P2 , P2 , (^2 
in terms of p2,2 and the previously determined functions. We obtain p2,2 with the 
integrability condition = 



(Mye) - Myr)) H,e-'^^'/^{l - i7,e"-=^/0 + W2 ^2^2] {ye) - [W2 [^^2] (l/) 



{Mve) - Myr)) (1 - i/ce-2-^^/-) - W2 [Qf,2] {ye) + [w^2 [Qf,2] {y) 



in which all functions are taken at argument 2iujc 
The component at frequency zero Q2,o obeys 



o = /:[g2,o] + P2,o + 



dy 2 dy"^ 



r ^ dy^ 2 dy^ r 



(75) 
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together with the boundary conditions 



Q2,o(ye, t) = 0, ^^(ye) = -P2,o(l -H)- 2H^ cosM/r) 



and a similar condition in t/r- 



Its general solution can be written 

Q2,o{y) = 



(ytoQo + Pto exp(-i/^) + P2fiQ-2"o{y) + Q2"oiy) y>yr 

(^2,oQo + p2,o exp(-?/2) + p2,oQ2"oiy) + Qt'M y < yr 



(76) 



where 



dy 4 dy"^ J 

and it is again convenient to construct the particular solution Q':^^ by differentiation 



QUy) 



1 — iujc dy 2{2 — iuc) dy"^ 
Gl d^Qo , H,G,{2 + u^^) d^Q, 



+ 



+ c.c. 



d^Qo 



l + ul dy^ (1 + cu2)(4 + cu2) dy^ 4(4 + ^2) dy^ 



(77) 



In this case, the four boundary conditions for (^2,0 are not independent and 
are not sufficient to determine the four unknowns ^2^05 '^2^05 /^2^05 /^2^o functions of 
lower order terms. This comes about because some choices of Q2,o are equivalent to 
changing the normalization of Qq. One should therefore eliminate them by imposing 
the condition /Z^(i?/Q2,o = 0. In this way, one obtains. 



P2,0 



ye 



1 

2UQ 



G^-Mf)ey'Jl^due-' 



ye 
yr 



(78) 



where 



7g(i/) = Gey + cos(u;c5/r) - uJcSm{uJc6/T) - 



Hc{2y' + 1; 



lH{y) = ^ycos{ujJ/T) - 2yUcSm{uJcS/T) + ^^^(^^ + ^) _ Hc{2y^ + 3y) 



-fi = -2{ye - yr)GcHc 



2 + u;! 



+ 



Hl{yl-yl) 



{l + ujl){A + ul) ' 4 + o;2 

(the notation [f]yl = fiye) — fiyr) is used). The derivatives of higher order of Q 
and (52,0; which are used in the following, can be obtained recursively by differenti 
ation of Eq. (|7|) and (^. 



2.2 
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B.2 Third order 



We can now proceed and study the third order terms. They obey the equation 



C[Q;]+n,it-6){G,^ + 







dt \ dy 2 dy 



,2 



2 



9?/ 2 9?/^ y Y (9?/ 2 ) 

together with boundary conditions 

Qz{ye) = 

^{ye) = -ns{t) + Hns{t-5)-2H^m{t-5)n,it-5) 
dy 

+H {ni{t)n2{t -6)+ni{t- 6)n2{t)) + if^n?(t -6)- H^ni{t)nl{t - 6) 

+ (H- H,)ni (t-6)- iJ5^e'"^(*-^)/" (80) 

dt 

and a similar condition holds at yr- 

The last two terms between brackets on the r.h.s. of (ITSp come from the antici- 
pation that it will be needed to have hi change on a slow time scale to cancel secular 
terms. The first term arises from the explicit time differentiation in Eq. ( ^OD and does 
not need special explanations. The second is less usual and comes from the delayed 
forcing //(t — 6) in (^). Formally introducing a slow time scale T = et, the delayed 
forcing is written i>{t — 6,T — e6). The second term between brackets in ([79| ) is pro- 
duced by the expansion to first-order in e — 5, T — e5) = u^t — S) — edri^it — 6) + ■ ■ ■. 
The last term in the boundary condition ( PU] ) appears in the same way. 

The forcing terms on the r.h.s. of (^) oscillate at frequencies ScUc and ujc- There- 
fore, we search Q3{y,t) and n^it) under the form 

Qsiy, t) = e^'-^'/^Q,M + e'^^'/^Q^M + c.c. 

n^{t) = e3^'^=*/"n3,3 + e^'^=*/"n3,i + c.c. (81) 

We focus on the terms at frequency Uc which are resonant with the first order terms. 
They obey the equation 

^ dy ^ 2 dy^ r dy + 2 dy^ j 
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\ dy 2 ay J 

The general solution of ( ^21) can be written 



Q^(y) = I ^sMy^ ^^c) + fit 02 (y, i^c) + ^3,l<0l + Qtl y>yr /g3^ 

^ I "3^01 (y, ^^^c) + PsMy^ i^c) + n3,iQi + Qi°i y <yr 



In the particular solution, Qi is the function that appears at first order, Eq. 
and as before, we can construct Q^^i by differentiation of lower order terms: 

Qtl = ^^^3,1 + ^1^3,1 + ^ll^ll'^3,l (84) 

where i is obtained from Qi by differentiation of 0i,2 and Qi with respect to A 

Ad ( ^ ^ i oifOxMy^^^c) + PtdxMy^^^c) + dxQliy,iujc) y > yr .^.x 
^'''^^^ I f3^dxMy, ^^c) + dxQl{y, lu,) y<yr, ^ > 

[ l + too, dy ^ 2{2 + tu;,) dy^ J' ^^^^ 

and 

Ac ^ iuJcS/r ( Gc dQ2,2 He d'^Q2,2 \ 

\l-iuJc dy ^2{2-iu:,) dy^ ) 

-ioJc5/r ( Gc dQ2fi He d'^Q2fi \ 



+ e 



1 + iujc dy 2(2 + iuc) dy"^ 



+ P'ifi\^c Qy + ^ Qy2 j+P2,2e [i + 2iuje dy ^4(l + za;J dy^ ^ 

Gc ^ Hed^\ _ , He fGed^Q, ^ H^ d'Q,\ 



l+uj^\ dy^ 3 dy^ J 4 + uj^\3 dy^ 8 dy^ 

^-2^..S/r Ge ( Ge d'Q{ H^ d'Qt ) 

l + iiUe\2{l + iLUe) dy"^ 2{3 + 2iuJe) dy^ ) 

2^..8|r He ( Ge d^Q\ ^ He d'Q\ \ 



e 



- p2,2e 



2{2 + iuJe) \3 + 2iuJe dy^ A{2 + iue) dy^ J 

2 + iUe ( Ge d'^Qo He d^Qo 



{\ — iuj (){\ ^ 2iijj e) \ 2^iLiJe dy"^ 2{3 + iuje) dy^ 



-iu^.^lrn '^+^^c f Ge d^Qo He d^Qo 

P2,oe ' Ge——. 7-— j-ir- + 



l+iuje \2 + iuje dy"^ 2{3 + iuje) dy'^ J 

^-^^^S/rjj 4 + / Ge dQl He d^Qp 

^^'^ '4(1 + iuOe){2 - iuoe) \3 + luJe dy^ 2(4 + iuje) dy^ ^ 
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He d'^Qo 



'A{2 + iuJe) yS + iUc dy^ 2(4 + icuj dy"^ J 

3 + iuJe ( Ge d^Qo He d^Q^' 



2(1 -iu;c)(l + \^3 + icu^ dy"^ 2{A + iuOc) dy^ 



+ e 



+ 



+ 



6 Vl+^c 4 + cj2 (l + icU,)(2 + ZCJe), 



A + iuJc dy'^ 2{5 + iuJc) dy^ 
6 + iujr 



Gc d'Qo^ He 



]{2 + iuOc)'^{2-iuJc) \h + iuOc dy^ 2{<o + iuOc) dy^ 



^7) 



Now, upon replacing 03 = one can try to determine a'^ , , [3'^ ,hs from the 
four boundary conditions on Q3,i{y). This provides a hnear inhomogeneous system 
for the four unknowns. The inhomogeneous terms are made from Q'i^iiy) and its 
derivatives evaluated at ye and yr- But there is a difficulty : since we are considering 
the resonant part of the third order terms, the linear operator coincides with the 4x4 
matrix obtained at first order which has been required to have a zero determinant. 
So, the equations for a^,P^,P^,h3 are solvable only if the inhomogeneous terms 
obey a solvability condition. In order to obtain it, we find it convenient to proceed 
as we did at linear order (see Eq. ( |63| , |65D). We obtain and in terms of h^ i and 



(53° from the 2x2 system given by the two boundary conditions at ye. We then obtain 
similar expressions for and (3^ — (3^ . Comparing the two obtained expressions 
for and requiring them to be identical provides the solvability condition. 



(88) 



where 



r dt 
^3 = -2iy2e--=^/-(p22 + p2o) + 3if;^e--^V- 



+ He 



P22(e-''"^'/" + e"^^"^) + p2o(l + e" 



-iuJc5/T\ 



~Hl{2 + e 



-2iujc^/T\ 



^9) 



With the help of Eqs. ( p^ , |85| , p6| , |87| ) , this gives the searched for equation of motion 
for ni 

= Aui- B\ui\^i>i (90) 



in which 



dT 



A 



-Wo 



Qli (ye) + W, a I (y) 



Mye) - Myr)) 



Qii] (ye) - [m [Qii] (y)]': + (Uye) - Myr)) '-He-^-^^/r 



Wo 



(91) 
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B 



Wo 



{ye) - \W2 \QI,\ {y)f_ + (Uve) - Uvr)) ^3 

Qii] (yo) - [m [Qii] {y)Y^ + (Mye) - Mvr)) ^ffe--^^/- 



(92) 



These expressions simplifies in the hmit 6/t ^ 0. In the particular case H = 0, one 
obtains Eq. (^3|) of the main text. 



C Effect of noise due to finite-size effects 



Inserting the noise in Eq. ([75|), we obtain 

dni 



dt 



Ahi - B\ni\^hi + D^Cif) 



(93) 



in which A and B are given by Eqs. (|9T1 , |92|) , while D is 



-W2 



D = r]- 



Qr 



iye) + 



Wo 



{y) 



Vr 



Wo 



Qii {ye 



where 



Qr 



W2 [Qii\ (y) 



(94) 



1 + iujc dy 

7] is given by Eq. (|27|) , and C is a complex white noise such that < C{t)C*{t') >= 
5{t - 1') 

The autocorrelation at zero time C(0) is given by 

C(0) = 1 + 2 < |ni(t)|2 > 

We deduce from Eq. ( pSf ) the Fokker-Planck equation describing the evolution of the 
p.d.f. of both real and imaginary parts of hi. This equation can be converted in an 
equation giving the stationary distribution Pr(p) of p = jiT-ip. It satisfies 



whose solution is 



IK 



Pr(p) 



d 



dp I dp 



— [[2Arp-2Brp' 



Ft 



exp 






/o°°exp( 


, |D|2-^ |D|2-^ , 


) dR 



and the autocorrelation at zero lag is 



C(0) 



J^Rexp {2j^R-^R^) dR 
/-exp(2^i?-^i?2)rfi? 



From this exact expression, it is not difficult to obtain the expressions (|29|j30| , |31|) of 
the main text. 
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From Eq. (PBD, one can compute the behavior of the autocorrelation function 
C{s). Far below the critical line, |ni| is small and the nonlinear term can be ne- 
glected. It is then easy to obtain Eq. (p2D of the main text. 

In the oscillatory regime far above the critical line, finite size effects provoke 
fluctuations of activity around the oscillation described by Eq. (|^). We consider a 
small perturbation, both in amplitude and in phase, of the 'pure' oscillation -hi 
ni(l + r) exp(z0). r is the perturbation in amplitude, while (j) is the perturbation in 
phase. To obtain the evolution equations for r and (p we apply standard stochastic 
calculus techniques (see e.g. Gardiner 1983, chapter 4), and obtain, 

rr = -Ar{2r + 3r^ + r^) + eCr + e \ ^ (95) 

2(1 + rj 

^^^'■.(2r + r2) + e-^ (96) 



Br " ' 1 + r 

in which e = \D\IR^ and C^j-i d are uncorrelated white noises. Note that the last term 
in the r.h.s. of Eq. (^) appears due to the fact that, upon discretizing Eq. ( P5| ) with 
a small time step dt, (f){t + dt) — (f){t) is of order y/dt, not dt. The calculation of the 
autocorrelation in terms of r and gives, keeping only the dominant term, 

C{s) = 1 + 2i?^ < cos {{uc + Auj)s/t + (pit + s)- (pit)) > . 

In order to calculate the autocorrelation we need to calculate the distribution of 
A0(s) = (p{t + s) — (pit). From Eqs. ( P3| , PUD we find that, to leading order in e, it has 
a Gaussian distribution with mean and variance 



l\s) 



\Dl_ 



s Bf r / 2ArS\ 2ArS 

- + 777^ 1 + 



_r 2B^Ar 

Averaging cos((u;c + Au)s/t + A0(s)) with such a distribution yields 
C(s) = 1 + 2i?2 cos ((u;c + Auj)s/t) exp (--f^{s)/2 



We find a damped cosine function as below the critical lines, but now the damping 
factor is no longer a simple exponential. For small times s -C r/ {ByB?), the damping 
is described by 

■ ^ ■ / \D?s\ 




exp ~ exp 




2 ) \ AR^T 

while for long times s ^ t/ [B^R^) 

( l\s)\ 
exp I — I ~ exp 

The damping time constant in both regimes is proportional to 1/|-Dp ~ N/C, i.e. to 
the inverse of the connection probability. When goes to infinity at C fixed the 
'coherence time' of the oscillation increases linearly with A^. 

The next order in e brings (after a rather tedious calculation) a small additional 
contribution to the variance, so that for long times 




exp — = exp 
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D Randomly distributed synaptic times 



The calculations performed in the case in which all synaptic times have the same 
value can be repeated in the more general situation in which synaptic times are drawn 
randomly and independently at each site with distribution Pr(5). The difference is 
that, in all equations were functions of 6 appears, we need to integrate these functions 
with the p.d.f.Pr(5). For example, we find that the critical line where the instability 
appears is given by 

(Myo) - MVr)) j Pr((5)e-"'^/^ci(5) = 

in which 

E Inhomogeneous networks 

We now relax the constraint that the number of connections received by a neuron be 
precisely equal to C. The connections are randomly and independently drawn at each 
possible site. They are present with probability C/N. In this situation, the dynamics 
of different neurons will depend on this number of connections they receive: this 
number is now a random variable with mean C and variance C(l — e). For example, 
their frequency will be a decreasing function of the number of connections. The 
connectivity matrix is defined by Jij = Jeij where for all i,j Cij = 1 with probability 
e. The distribution of frequencies in the stationary state in such a situation has been 
obtained, for the case of a network with both excitatory and inhibitory neurons, by 
(Amit and Brunei 1997b). The distribution of stationary frequencies can be obtained 
as a special case of this analysis. We briefly recall here the main steps of this analysis, 
before turning to the stability analysis. 

Averaging the synaptic input only on the randomness of spike emission times of 
presynaptic neurons, we get that the mean and the variance of local inputs are given 

by 

j 3 

Since the number of inputs to each neuron is very large, the spatial distribution 
of the variable J2j^ij^j, which determines completely the spatial distribution of 
and (7, will be close to a Gaussian whose two first moments can be calculated as a 
function of the two first moments of the spatial distribution of frequencies: 

j 



Ql {ye 



Wo 



Q\ {y) - (97) 
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Thus the variable _ 

J2j eijUj - Cu 



Zi 



— ez/^ 



has a Gaussian distribution, p{z) = exp{—z^/2)/\^2n. Thus a neuron receives, with 
probabihty p{z), a local input with moments 



p,{z) = -Jt{Cv + z^jc (z/2 - ev^) ) (99) 

and 



a\z) = rT{Cv+zJC[u^-tV^)) (100) 



E.l Distribution of frequencies in stationary state 

In the stationary state the frequency of a neuron with moments p{z) and (j{z) is 
given by 

Uo{z) = [r^ duexp{u'){l + erf(«)) j (101) 

The two first moments of the distribution of frequencies can then be determined in 
a self-consistent way, using 

Vo= dzp{z)uo{z), ul= dzp{z)ul{z) 



These equations, together with Eqs. (|99|J100| , p!oTD , fully determine the whole distri 



bution of stationary frequencies, which can be obtained using the relation 

P{v) = I dzp{z)6{u - uo{z)) 



E.2 Linear stability analysis 

The linear stability analysis of Section ^ can be generalized to the inhomogeneous 
network. We give here the main steps of this analysis. 

We expand the frequencies around the stationary frequency, 

u{z) = Uo{z) (1 + ni{z, t) + ...), 

and, defining for each z y = {x — po{z)) / ao{z) , 

p = ^^^my,^) + Qiiy,^,t) + ---) 

The moments of the spatial distribution of frequencies can be expanded in the same 
way, 

17 = z7o (1 +ni(t) + ...), 
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where 



ni{t) = — dzp{z)iyQ{z)ni{z,t) 
z/q J 

nlit) = ^ dzp{z)ul{z)ni{z,t) 
The Fokker-Planck equation at first order is 



r- 



dt 



H^{z)n,{t -5) + H,{z)nl{t - 6)) d^Q^ 



+ (Gi{z)n^{t -6)+ G2{z)nl{t - 6) 



dy 



where 



JCuqt — eJrzJ C {vq — ez/p) = 



cro(^) 



Goiz) 



JrzjC[ .l-evt)^^ 
2ao(z) 



Hi{z) 



J-Cv,r-e.PrzJC[ul-eVl)=J^.^ 



(TniZ 



H,(z) = 



2aiiz) 



The eigenmodes of Eq. (|102|) can be written 

Qi{y, z,t) = Qi{y, z) exp^iut/r) + c.c. 

ni{z, t) = ni{z) exp{iut/T) + c.c. 
leading to the solvabihty conditions, for each z 

ni{z) = I{z)hi + J{z)hi 

where 



(102) 



W2 [Ri] (ye) - W2 [Ri] {y) _ + H,{z)e-''^'l^ Mye) - HVr) 

I{z) = ^ = ^ : 

Mve) - MVr) 

W2 [R2] (ye) - \W2 [R2] {y)Y^ + H2{z)e-'^'/- (Mve) - MVr) 
J{z) = ^ = ^ 



with 



Rio = e 



( GiM dQojy) ^ Hi^2{z) d^Qojyy 
\l + iw dy 2{2 + iw) dy'^ 



(103) 
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Multiplying the above equation by puo (2pi/q) and integrating with respect to z we 
obtain 

— < VqI >— < Z/qJ >^ 

rii — — — — rii H — — nf 

hi = 2 ° ni + 2 ° hi 

where we use the notation < . . . >= / dzp{z) .... The instability point together with 
the associated frequency are given by the condition that the associated determinant 
vanishes, i.e. 

_ < vqI > ^ < vlJ > ^ < >< VqJ > - < vp >< VqI > 
J-— — +^ — ^ +^ — 

The relative degree of synchrony of population z with the collective oscillation is 
given by 

h,{z) = h; i{z) + j{z)- — ^ 
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